Exploration of real bulk RNA-seq datasets

GFRN signature dataset

Introduction

The growth factor receptor network (GFRN) contributes to breast cancer progression and drug response. This RNA-Seq dataset (Rahman et al. 2017) is designed to develop gene signatures for several GFRN pathways: AKT, BAD, HER2, IGF1R, RAF1, KRAS, and EGFR. The study used recombinant adenoviruses to express these genes in case samples and produce green fluorescent protein (GFP) in control samples, using replicates of human mammary epithelial cells (HMECs). RNA-Seq data are collected from these HMECs overexpressing GFRN genes and GFP controls.

Overall design

The dataset contains 3 batches, with a total of 89 samples:

• Batch 1 (GSE83083): 41 samples. 6 replicate samples of each for AKT, BAD, IGF1R, and RAF1, 5 replicates for HER2, and 12 replicates for GFP controls
• Batch 2 (GSE59765): 12 samples. 6 replicates of each for EGFR and its corresponding controls
  
• Batch 3 (GSE83083): 36 samples. 9 replicates of each for 3 types of KRAS mutants and GFP control

Load data & basic info

A few basic statistics of the dataset:

######  Samples  #######
- Total number of samples:
[1] 89

- Number of samples in condition groups:
group_num
 0  1  2  3  4  5  6  7  8  9 
27  6  6  6  5  6  9  9  9  6 
In batch 1:
  akt   bad gfp18  her2 igf1r   raf 
    6     6    12     5     6     6 
In batch 2:
        egfr gfp_for_egfr 
           6            6 
In batch 3:
 gfp30 krasgv krasqh kraswt 
     9      9      9      9 

- Number of samples in batches:
batch
 1  2  3 
41 12 36 


######  Genes  #######
- Total number of genes:
[1] 23368

- Number of genes with 0 counts:
10097 43.2%

- Number of genes with maximum 10 counts:
7868 33.7%

- Number of genes with ONLY 0s:
2500 10.7%


######  Counts  #######
- Range of count in whole dataset:
[1]       0 2270900

- Average count in whole dataset:
[1] "Mean: 886.527"
[1] "Median: 45"

For the following analysis, we remove genes with only zero values across all samples.

## Remove genes with all 0 counts
sigdata <- sigdata[apply(cts_mat, 1, function(x){!all(x==0)}), ]
cts_mat <- assay(sigdata, "counts")
rownames(cts_mat) <- paste0("gene", 1:nrow(sigdata))
print(dim(cts_mat))

[1] 20868    89

The figure below shows distribution of percentages of zero counts for the remaining genes (each data point represents a gene, i.e., how many 0 counts are there in each gene).

######  Proportion of genes with X% zeros across samples  #######
 <10.0%  <20.0%  <30.0%  <40.0%  <50.0%  <60.0%  <70.0%  <80.0%  <90.0% <100.0% 
  0.703   0.045   0.020   0.016   0.017   0.017   0.020   0.025   0.039   0.097 

Existence of batch effect

We first visualize the data using PCA and hierarchical clustering, to see the batch effect in the dataset. Especially, Batch 2 separates from the other two batches in PCA plots.

PCA

Dendrograms

Library size

Next, we take a look at the observed library sizes of samples, and to see if there is a difference in library size across batches.

Range of observed library sizes across samples:
[1] "1.0e+07" "3.4e+07"

Gene expression distribution

Our next step is to interrogate the gene expression distributions.

We first simply plot the counts, CPM, and log2CPM of 20 randomly chosen genes in the dataset. These figures are stored under directory specified by output_dir.

Fold changes in mean

We then use edgeR to perform differential expression analysis, and calculate (log-) fold changes of all genes between condition groups and batches.

Change WRT condition

To rule out the impact of batch effect, we limit our DE analysis with regard to biological condition in one of the batches. Here we specify the batch to be Batch 2: EGFR vs GFP control.

####  Differential genes relative to biological condition 
####  (within batch specified by 'test_batch')
test_batch <- 2
counts_batch <- cts_mat[, batch==test_batch]
cond_batch <- group_num[batch==test_batch]
y <- DGEList(counts=counts_batch, group=as.factor(cond_batch))
y <- calcNormFactors(y, method="TMM")
design <- model.matrix(~as.factor(cond_batch))
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, coef=2)
de_res <- topTags(qlf, n=nrow(sigdata))
 
up_genes <- rownames(de_res)[which(de_res$table$logFC >= 0)]  # all up-regulated genes (descending order in terms of p value)
head_up_genes <- head(up_genes, n=50)  # top-50 significant up-regulated genes
tail_up_genes <- tail(up_genes, n=50)  # 50 least significant up-regulated genes
cond_up_inds <- up_genes 
down_genes <- rownames(de_res)[which(de_res$table$logFC < 0)]  # all down-regulated genes
head_down_genes <- head(down_genes, n=50)  # top-50 significant down-regulated genes
tail_down_genes <- tail(down_genes, n=50)  # 50 least significant down-regulated genes
cond_down_inds <- down_genes

######  Changes in terms of average gene counts  ######
     top50 up top50 down 100 least DE
GFP    1535.8     4700.0        325.2
EGFR   2142.5     3570.9        321.2

######  Range of fold change across biological conditions  ######
          top50 up top50 down 100 least DE
min FC       1.421      0.296        0.988
median FC    1.652      0.648        1.000
mean FC      1.809      0.616        0.999
max FC       3.657      0.766        1.000

######  Range of Log(fold change) across biological conditions  ######
             top50 up top50 down 100 least DE
min logFC       0.352     -1.219       -0.012
median logFC    0.502     -0.434        0.000
mean logFC      0.593     -0.484       -0.001
max logFC       1.297     -0.266        0.000

Conclusions:

• Biological fold change comparing average counts in top DE genes is about 1.3-1.4.
• Gene-wise fold change ranges in top DE genes from 1.7-1.8.

Change WRT batch

We are also interested to see how the average (mean) gene counts change across batches (i.e. degree of mean batch effect).

To explore this mean batch effect, we first perform DE with regard to the batch variable. To rule out the impact of biological condition, we perform this analysis within control samples only. Since from the PCA plot, Batch 2 separates from the other two batches more, we compare Batch 1 and 2 in this case.

####  Differential genes relative to batch 
####  (within condition 0 / condition specified by 'test_cond')
test_cond <- 0
counts_cond <- cts_mat[, (group_num==test_cond) & (batch!=3)]
batch_cond <- batch[(group_num==test_cond) & (batch!=3)]
y <- DGEList(counts=counts_cond, group=as.factor(batch_cond))
y <- calcNormFactors(y, method="TMM")
design <- model.matrix(~as.factor(batch_cond))
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, coef=2)
de_res <- topTags(qlf, n=nrow(sigdata))
up_genes <- rownames(de_res)[which(de_res$table$logFC >= 0)]  # all up-regulated genes
head_up_genes <- head(up_genes, n=50)  # top-50 up-regulated genes
tail_up_genes <- tail(up_genes, n=50)  # 50 least up-regulated genes
batch_up_inds <- up_genes
down_genes <- rownames(de_res)[which(de_res$table$logFC < 0)]  # all down-regulated genes
head_down_genes <- head(down_genes, n=50)  # top-50 down-regulated genes
tail_down_genes <- tail(down_genes, n=50)  # 50 least down-regulated genes
batch_down_inds <- down_genes

######  Batch differences in terms of mean gene counts  ######
       top50 up top50 down 100 least changed
Batch1    389.6     4653.5              42.6
Batch2   4001.1      158.2              57.3

######  Range of fold change between batches  ######
- All genes - FC:
0 34381.28

- All genes - log(FC):
-15.972 10.445

- Top genes - FC:
           top50 up top50 down 100 least DE
min FC        7.664      0.000        0.979
median FC    27.103      0.001        1.000
mean FC     404.109      0.008        0.996
max FC    12757.952      0.059        1.000

- Top genes - log(FC):
             top50 up top50 down 100 least DE
min logFC       2.037    -15.972       -0.021
median logFC    3.300     -6.866        0.000
mean logFC      6.002     -4.881       -0.004
max logFC       9.454     -2.836        0.000

Conclusion:

• Mean batch effect fold change comparing average counts in top genes is quite large. But remember, these are top genes that change across batches.

Batch mean

Median (of all genes) gene-wise mean:

         Batch1 Batch2 Batch3
Controls  90.17  98.33 114.94
Cases     80.33  95.67 151.22

Scaled by that of batch 1:

         Batch1 Batch2 Batch3
Controls      1   1.09   1.27
Cases         1   1.19   1.88

Batch dispersion

Aside from average gene counts, we are interested in comparing dispersions in the batches. Difference in the dispersions across batches suggests variance batch effect that are not accounted for by Poisson model. We visualize the gene-wise dispersion estimates from edgeR:

####  Statistics about dispersion estimates from the batches  ####
             Batch1 Batch2 Batch3
min disp.    0.0001 0.0001 0.0001
median disp. 0.0075 0.0059 0.0027
mean disp.   0.0235 0.0309 0.0147
max disp.    1.0377 2.8392 1.6146

Conclusions:

• Max mean dispersion difference: Batch 2 is about 2 times of those in Batch 3.

Batch variance

Median (of all genes) gene-wise variance:

         Batch1 Batch2 Batch3
Controls 488.66 475.07 518.43
Cases    179.53 870.22 671.43

Scaled by that of batch 1:

         Batch1 Batch2 Batch3
Controls      1   0.97   1.06
Cases         1   4.85   3.74

Impact of 0 counts on gene distributions

We have seen that a lot of genes in the dataset contain some portion of zeros. Having 0 counts can have a great impact on mean and dispersion estimation. We evaluate such impact in this section, and see how NB distribution fits to non-zero portion of genes. We first take the non-zero portion of the genes.

Fold change in mean in non-zero portion

Change WRT condition in batch 2

Within batch 2, we compare mean (here refers to non-zero mean) across biological condition groups. Specifically, we saved the previous DE analysis results using edgeR, on counts containing zeros. We call all genes with \(\log(FC) >= 0\) as “up” genes, and \(\log(FC) < 0\) as “down” genes. The top 50 “up” and “down” genes (with smallest p values) are still determined by the previous DE containing zeros.

######  Changes in MEAN non-zero gene counts  ######
     top50 up top50 down 100 least DE
GFP    1535.8     4700.0        668.3
EGFR   2142.5     3570.9        673.9

######  Changes in MEDIAN non-zero gene counts  ######
     top50 up top50 down 100 least DE
GFP    1078.5     1691.5        178.5
EGFR   1499.0     1188.5        175.5

Conclusions:

• Biological fold change (in non-zero portion) comparing average counts in top genes is about 1.3-1.4 (mean), 1.39 - 1.42 (median).**
• O counts did not strongly affect biological signal in terms of average counts in top genes.

Change WRT batch in control samples

Similarly, we compare non-zero mean in top genes across batches within control samples. Here “top genes” refer to top gene lists from previous DE within control samples comparing batches.

######  Changes in MEAN non-zero gene counts  ######
        top50 up top50 down 100 least DE
Batch 1    394.8     4653.5        184.7
Batch 2   4001.1      189.9        213.5

######  Changes in MEDIAN non-zero gene counts  ######
        top50 up top50 down 100 least DE
Batch 1    176.5       2048            2
Batch 2   2720.5         41            2

Conclusions:

• Mean batch effect fold change comparing average counts in top genes is very large.
• O counts did not strongly affect mean batch effect in terms of average counts in top genes.

Batch mean with 0 counts removed

Median (of all genes) gene-wise mean after removing zeros:

         Batch1 Batch2 Batch3
Controls 138.96 217.17 184.39
Cases    158.58 214.17 223.83

Scaled by that of batch 1:

         Batch1 Batch2 Batch3
Controls      1   1.56   1.33
Cases         1   1.35   1.41

Batch dispersion with 0 counts removed

Dispersion estimates are more likely impacted by 0 counts. We re-estimate the dispersions in each batch, using non-zero portion of each gene. Then we compare the dispersion estimates across batches.

####  Statistics about dispersion estimates (from NON-ZERO portions) from the two batches  ####
             Batch1 Batch2 Batch3
min disp.    0.0000 0.0000 0.0000
median disp. 0.0268 0.0512 0.0148
mean disp.   0.0271 0.0522 0.0163
max disp.    0.8216 2.2998 0.6972

Conclusions:

• The absolute value of maximum dispersion are reduced after 0 counts are removed.
• In this report, dispersion are calculated by different methods before and after removing 0 counts. Not directly comparable.
• We show in another analysis (gfrnSig_run_nonzero_disps.R) that the relationship between batches did not change much after removing 0s.

Batch variance with 0 counts removed

Median gene-wise variance after removing zeros:

          Batch1  Batch2  Batch3
Controls 1475.48 3331.57 1758.75
Cases     657.83 7140.48 2179.06

Scaled by that of batch 1:

         Batch1 Batch2 Batch3
Controls      1   2.26   1.19
Cases         1  10.85   3.31

Is there a difference in zero portions across batches?

Another question related to batch effect, is whether there is any difference in the percentage of 0 counts in a gene between batches.

Conclusions:

• There is a certain level of difference in zero fractions between batches. Specifically, Batch 3 has more relatively “complete” (fewer 0 counts) genes than Batches 1 & 2.

Goodness of fit for NB distribution to non-zero portion of genes

It is difficult to evaluate quantitatively the goodness-of-fit of NB distribution on data. Instead, we generate Q-Q plots for 100 randomly selected genes, comparing data with NB distribution using ML estimated parameters. We generate these plots within batch 2, for cases and control samples separately. Figures are stored under {output_dir}/qqplots/.

More justification for having different dispersion parameter for different batches

Finally, a very important assumption of ComBat-seq model (compared to adding batch covariates in linear model of edgeR and DESeq), is that having a different dispersion parameter for each batch is beneficial for DE performance. As argued by many, if not all, of these papers, accurate estimation of dispersion paramters are important for DE detection. As we have seen above, there is some difference in dispersion distribution across batches. If dispersion estimates within each batch and using the whole study are not consistent with each other, then it is justification that we may need different dipsersion parameter for different batches, and it could be good for DE.

Here we generate scatter plots comparing dispersion estimates within each batch to estimates from whole dataset (including batch as a covariate in whole study).

Gene-wise dispersion estimates

Trended dispersion estimates (log10 scale)

Session info

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] MASS_7.3-51                 ggdendro_0.1-20             scales_1.0.0                gridExtra_2.3              
 [5] reshape2_1.4.3              ggplot2_3.0.0               dendextend_1.9.0            DESeq2_1.20.0              
 [9] SummarizedExperiment_1.10.1 DelayedArray_0.6.6          BiocParallel_1.14.2         matrixStats_0.54.0         
[13] Biobase_2.40.0              GenomicRanges_1.32.7        GenomeInfoDb_1.16.0         IRanges_2.14.12            
[17] S4Vectors_0.18.3            BiocGenerics_0.26.0         edgeR_3.22.5                limma_3.36.5               

loaded via a namespace (and not attached):
 [1] bitops_1.0-6           bit64_0.9-7            RColorBrewer_1.1-2     prabclus_2.2-6         tools_3.5.1           
 [6] backports_1.1.2        R6_2.3.0               rpart_4.1-13           Hmisc_4.1-1            DBI_1.0.0             
[11] lazyeval_0.2.1         colorspace_1.3-2       trimcluster_0.1-2.1    nnet_7.3-12            withr_2.1.2           
[16] tidyselect_0.2.5       bit_1.1-14             compiler_3.5.1         htmlTable_1.12         labeling_0.3          
[21] diptest_0.75-7         checkmate_1.8.5        DEoptimR_1.0-8         mvtnorm_1.0-8          robustbase_0.93-3     
[26] genefilter_1.62.0      stringr_1.3.1          digest_0.6.18          foreign_0.8-71         XVector_0.20.0        
[31] base64enc_0.1-3        pkgconfig_2.0.2        htmltools_0.3.6        htmlwidgets_1.3        rlang_0.3.0           
[36] rstudioapi_0.8         RSQLite_2.1.1          bindr_0.1.1            mclust_5.4.1           acepack_1.4.1         
[41] dplyr_0.7.7            RCurl_1.95-4.11        magrittr_1.5           modeltools_0.2-22      GenomeInfoDbData_1.1.0
[46] Formula_1.2-3          Matrix_1.2-14          Rcpp_0.12.19           munsell_0.5.0          viridis_0.5.1         
[51] stringi_1.2.4          whisker_0.3-2          yaml_2.2.0             zlibbioc_1.26.0        flexmix_2.3-14        
[56] plyr_1.8.4             grid_3.5.1             blob_1.1.1             crayon_1.3.4           lattice_0.20-35       
[61] splines_3.5.1          annotate_1.58.0        locfit_1.5-9.1         knitr_1.20             pillar_1.3.0          
[66] fpc_2.1-11.1           geneplotter_1.58.0     XML_3.98-1.16          glue_1.3.0             latticeExtra_0.6-28   
[71] data.table_1.11.8      gtable_0.2.0           purrr_0.2.5            kernlab_0.9-27         assertthat_0.2.0      
[76] xtable_1.8-3           viridisLite_0.3.0      class_7.3-14           survival_2.42-6        tibble_1.4.2          
[81] AnnotationDbi_1.42.1   memoise_1.1.0          bindrcpp_0.2.2         cluster_2.0.7-1       

References

Rahman, Mumtahena, Shelley M MacNeil, David F Jenkins, Gajendra Shrestha, Sydney R Wyatt, Jasmine A McQuerry, Stephen R Piccolo, et al. 2017. “Activity of Distinct Growth Factor Receptor Network Components in Breast Tumors Uncovers Two Biologically Relevant Subtypes.” Genome Medicine 9 (1). BioMed Central: 40.